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We study the Landau damping of ferromagnetic magnons in Fe, Co, and Ni as the dimensionaUty 
of the system is reduced from three to two. We resort to the ab initio hnear response time dependent 
density functional theory in the adiabatic local spin density approximation. The numerical scheme 
is based on the Korringa-Kohn-Rostoker Green's function method. The key points of the theoretical 
approach and the implementation are discussed. We investigate the transition metals in three 
different forms: bulk phases, free-standing thin films and thin films supported on a nonmagnetic 
substrate. We demonstrate that the dimensionality trends in Fe and Ni are opposite: in Fe the 
transition from bulk bcc crystal to Fe/Cu(100) film reduces the damping whereas in Ni/Cu(100) 
film the attenuation increases compared to bulk fee Ni. In Co, the strength of the damping depends 
relatively weakly on the sample dimensionality. We explain the difference in the trends on the basis 
of the underlying electronic structure. The infiuence of the substrate on the spin-wave damping 
is analyzed by employing Landau maps representing wave-vector resolved spectral density of the 
Stoner excitations. 
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I. INTRODUCTION 

The properties of excited magnetic states are of great importance in the fundamental and applied magnetism. 
Their spectrum determines the thermodynamical properties of magnets, including the Curie temperature^"— . The 
excitations contribute to the electronic specific heat^ and the electrical and thermal conductivity, couple to charge 
degrees of freedom^"—, limiting the life times and the mean free paths of excited electrons^"— , and can provide a 
coupling mechanism in high temperature superconductors alternative to phonons^Sri^^. The control of spin dynamics 
and its attenuation are the central problems in the rapidly growing field of spintronics^"— . 

Magnetically ordered materials feature a class of collective spin-flip excitations called spin-waves or magnons^—. 
We focus on them in this paper. The spin-waves are usually pictured as a coherent precession of atomic moments 
around the direction of ground state magnetization'^^ . Every magnon carries certain energy ujq and crystal momentum 
q and changes the magnetization of the system by 2/xb . Metallic magnets feature another type of spin- flip excitation 
termed Stoner excitations. They are electron transitions between one-electron states with opposite spin projections. 
The precessing magnetization can flip spins of single electrons, creating Stoner pairs, which leads to the damping of the 
moments' precession. The attenuation of collective excitations due to the interaction with single-particle continuum 
was first considered by Landau^i^^ and is commonly referred to as Landau damping. This mechanism of attenuation 
is dominant in metals that are of our primary interest here. Other decay processes, which could be captured in a 
language of magnon-magnon, magnon-phonon or magnon-electron interaction, to mention the few, are not considered 
in this paper. 

The theoretical description of spin dynamics that accurately treats both spin-wave and Stoner excitations has 
proven to be a challenge. Among different formalisms, a particularly transparent approach is offered by the linear 
response time dependent density functional theory^^ , relying on the evaluation of wave- vector dependent dynamic spin 
susceptibility. Such calculations are, however, very demanding both from the point of view of algorithmic complexity 
and computer resources and for a long time they were restricted to simple bulk systems. "^^'^^ Recently, we have 
developed an efficient numerical scheme allowing to evaluate the spin susceptibility of complex magnets and applied 
it to study energies and life-times of magnons in complex bulk phases^ and ultrathin films^i^. In Sec. Ill A] we give 
a brief historical overview of the developments in this field and relate our calculational approach to works of other 
authors. 

The properties of three transition metal ferromagnets Fe, Co and Ni have been for many years a subject of intensive 
theoretical studies (see, e.g.^^). Much attention has been also devoted to understanding spin- waves in the ultrathin 
films absorbed on metallic non-magnetic substrates^i^, as such systems are of enormous practical importance. In the 
pioneering theoretical works^"— it was argued that the Landau damping of magnons in the ultrathin films should be 
generally more severe than in the corresponding bulk.*^ The broken out-of-film-plane translational symmetry should 
increase the cross section for scattering of magnons and Stoner excitations. Additionally, the non-magnetic substrates 
feature infinite continuum of gapless Stoner excitations. Muniz et al.'^^ pointed to the Stoner transitions involving 
the electronic states of the substrate as an important contribution to the Landau damping. We have shown^^ that 
such simple dimensionality arguments are insufficient and the actual magnitude of Landau damping is sensitive to the 
details of the hybridization between the electron states of the film and the substrate. 

In this paper, our goal is two-fold. First, we provide a thorough expostion of our theoretical approach and outline 
key details of its numerical implementation. Second, we report detailed study of the influence of the dimensionality 
of the system on the spin-wave damping in iron, cobalt and nickel. We show that the trend in the damping variation 
upon the transition from bulk to film is opposite for Fe and Ni whereas in Co the attenuation is practically unchanged. 
We provide an explanation of these trends on the basis of the properties of the electronic structure. 

The paper is organized as follows. In Sec. |TT]we delineate the formalism of the time dependent density functional 
theory for spin-flip excitations in the linear regime. Sec. IIIII covers the details of the computer implementation. The 
results of the numerical studies are gathered and discussed in Sec. IIVI 
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II. FORMALISM 
A. Historical perspective 

Until recently, the main body of the theoretical studies of spin-waves was based on the adiabatic treatment of 
magnetic degrees of freedo m^ where one maps the spin system onto an effective Heisenberg Hamiltonian of atomic 
moments which is much easier to study. The approach has acquired many different forms in the literature, namely 
the magnetic force theorem (MFT)^, frozen magnon techniqu o^^i^^'^^ and static transverse susceptibility method 
(STSM)^S. Bruno^ has shown that the MFT machinery in the form originally suggested gives a systematic error in 
the estimated Heisenberg exchange parameters. STSM method based on a consistent account for the external field is 
free from this error. 

These adiabatic methods utilize density functional theory and, therefore, do not involve adjustable parameters. 
They are exact in ferromagnets in the limit of small magnon momenta, i.e. they allow to determine correctly the 
spin-wave stiffness constant. In this approach, however, the presence of the Stoner excitations is neglected and no 
prediction regarding magnon life-times can be made. Furthermore, for higher energies one expects that the coupling 
to the Stoner continuum, apart from leading to the decay of spin-waves, can also cause a renormalization of the 
magnon energy.^^ In the latter case the adiabatic description must fail also in the prediction of the dispersion. 

As mentioned in the Introduction, the damping is captured in the calculations of wave-vector and frequency- 
dependent transverse magnetic susceptibility x(Q: ^) where spin- waves and Stoner excitations are treated on an equal 
footing. One associates the spin-wave excitations with the poles of x(q, w) in the complex energy plane. Roughly, 
the real part of the pole position gives the energy of the spin-wave and the imaginary part can be identified with the 
inverse life-time of the excitation. 

Despite the commonly accepted importance of the dynamic susceptibility approach, very few ab initio calculations 
along this line exist due to the complexity of corresponding numerical algorithms and high numerical costs of such 
calculations. The calculations performed so far deal mostly with elementary bulk systems. In 1980 Cooke et a/.— 
computed the susceptibility of Ni and Fe using the random phase approximation to the susceptibility starting from 
an idealized band structure. The dynamic method became particularly powerful after the parameter free linear 
response time dependent density functional theory (LRTDDFT) was developed in 1985 by Gross and Kohn,??. (An 
earlier work of Callaway and Wang^ contains already all key ingredients of the magnetic response in the local density 
approximation.) From this time on, several works addressing the paramagnetic susceptibility appeare d'^^i^^" — taking 
advantage of the formalism. In 1998 Savrasov^ presented the first calculation of the susceptibility using linear 
response density functional theory for magnetically ordered systems. During the past few years we developed a 
novel and efficient computer implementation based on the Korringa-Kohn-Rostoker (KKR) Green's function (GF) 
method^Sr— to address the spin susceptibilities of complex bulk systems and thin magnetic films. 

Recently, Lounis et alM-t^ reported on a considerable development of the earlier used empirical tight-binding scheme 
for the calculations of the transverse dynamic magnetic susceptibility, see e.g4^ that takes advantage of KKR Green's 
functions. Several simplifications were introduced with respect to the standard KKR method such as the neglect of 
the energy dependence of the electronic wave functions as well as the restriction to the d-electron states only. 

In our studies we rely on the treatment of the Kohn-Sham Hamiltonian of the DFT theory without simplifications 
beyond those used in the standard DFT calculations. In order to represent the response functions we use complete and 
numerically efficient Y-Ch basis^ consisting from the products of the spherical harmonics and Chebyshev polynomials. 
In particular, we do not need to assume the rigid rotation of atomic magnetic moments^ in the description of the 
spin response. Our experience shows that the properties of the spin- wave damping considered here depend sensitively 
on the structure of the Stoner continuum and therefore on the detailed form of the wave functions. 

It is important to mention a number of works on the calculation of the dynamic magnetic susceptibility performed 
within the framework of the many body perturbation theory (MBPT) This approach has a strong potential for the 
study of spin-wave excitations, especially in the systems with strong electron correlations where MBPT approach can 
be advantageous with respect to the DFT theory in the local density framework. However, up to now, the MBPT was 
applied only to elemental bulk ferromagnets^^'^'' and time will show if the computational complexity of the scheme 
will allow the efficient study of nanoscopic systems. It should be also noted that the practical implementations 
of the method necessarily simplify original Hedin's equations, e.g., by retaining only the ladder diagrams in the 
T-matrijsii^— i2£ used to construct spin-spin correlation function. The energy dependence and non-locality of the 
screened Coulomb interaction is often neglected as well. At this level of approximations the MBPT yields essentially 
the same physical picture of magnons and their damping as the adiabatic local spin density approximation (ALSDA) 
commonly used in time dependent density functional theory (TDDFT). 

Concerning earlier theoretical studies of the spin- wave damping in the two-dimensional magnetic systems, up to our 
knowledge, there is only one series of investigations based on an empirical tight binding scheme and random phase 
approximation to the dynamic susceptibility These pioneering works yielded rich qualitative information on 
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the spin-flip dynamics as well as on the relation between the parameters of the model and physical observables. 
Unfortunately, the strong sensitivity of the results to the parametrization of the underlying electronic stricture limits 
the usefulness of this approach. 

In the context of comparing different calculational techniques it is important to address the question of preserving 
the spin rotational invariance. In the absence of magnetic anisotropy the rigid rotation of spin system costs no 
energy. In magnetically ordered materials this implies that the spin-wave spectrum becomes gapless and in the limit 
of the vanishing momentum the magnon energy should tend to zero, featuring a Goldstone modei ''^i^^ Formally, both 
LRTDDFT and MBPT approaches satisfy the conditions of the Goldstone theorem, albeit the concrete numerical 
realizations lead to deviations of the magnon energy from zero in the limit oi q — 0. A method to compensate for 
numerical errors and to recover the Goldstone mode was recently discussed in Refsj^i^. Taking into account the 
importance of the issue and the absence of full agreement with the conclusions Lounis et al. arrived at, we address 
the problem in some details in Sec. IIIDI 

B. Linear response time dependent density functional theory 

The theoretical description of the time evolution of magnetization poses a complex problem in the solid state 
physics. Time dependent density functional theory^^ offers a tool for a parameter free description of time evolving 
system, also under the influence of dynamic external flelds. The method is capable of describing arbitrarily large 
excitations, including non-linear effects like higher harmonic generation, but at present it is computationally too 
expensive to be efficiently applied to solids. In 1985 Gross and Kohn laid foundations of linear response time dependent 
density functional theory framework'^^, which is applicable when the time dependent perturbation is small. The latter 
formalism allows to determine directly the density-density response function starting from the knowledge of Kolm- 
Sham (KS) ground state eigensystem. The knowledge of the response function is sufficient to characterize the excited 
states of the spin system. It is well suited for the ab initio description of magnetic excitations in periodic solids. 
Our efforts in this field concentrated on the development of an accurate and effective calculational scheme applicable 
to complex magnets with many inequivalent atoms. In this and the next sections we expose the formalism and our 
numerical approach in some details emphasizing the aspects we consider important for the accuracy and efficiency. 
We resort to the adiabatic local spin density approximation ( ALSDA) j36^58,75- _77 unless specified otherwise Rydberg 
atomic units are used throughout. 

We make use of the field operators^ 

= ^0j(xq)cj, ^^(xa) = ^0j(xa)*c], (1) 

where xa denote spatial and spin degrees of freedom, stands for a complete orthonormal set of single particle 
orbitals, and fermionic operators Cj and cj, respectively, destroy and create a particle in the orbital j. The particle 
and magnetization density operators read 

aXx)^<^$t(x^)^(x/3), (2) 

where cr^^ = Sap corresponds to the density operator and cr^'^'^ = a^^v^^ are standard Pauli matrices. Einstein 
summation convention is used for a, /3, . . . indices. 

We consider a time dependent perturbation Hamiltonian including oscillating electric potential and magnetic field 
coupled exclusively to the spin degree of freedom 

3 

H'^'^it) = J2 dxa*(xt)S,(xt). (3) 

i=a •' 

Relativistic effects and diamagnetic response are neglected. All many-body operators in this work which contain time 
variable are assumed to be in the Heisenberg picture. The driving potential is defined as a four-vector 

5(xt) = (-|e|V(xt),-^B(xO), (4) 

where V stands for the external electric potential, e is the charge of electron, B is the external magnetic field, ^-q 
denotes Bohr magneton and g is the electron g-factor. The approximation g — 2 \s used in this work. We define the 
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charge and magnetization density induced by the external time-dependent potential as 

5n\Kt)^{a\^))Jt)~n\^), (5) 

where n*(x) = ((t'(x)) is the unperturbed ground state charge and magnetization density and {o)^^{t) stands for 
the expectation value of operator o at time t when the perturbation is active. 6n^ is related to Sj(xt) through the 
retarded density-density response function 

x^^(x,x', t-t') = -x0it - t'){[a\^t),a^{^'t')]), (6) 

where [A^ B] = AB — BA. In the frequency domain one obtains 

3 

feXx,w) = j dU''^'Sn\x,t) = J dx'x^^(x,x',cj)Sj(x',cj). (7) 

LRTDDFT allows to compute the generalized susceptibility in the following two step procedure<^i2S First, one 
considers the Kohn-Sham susceptibility 

f ' \ ^ j ff f , (/)fc(xa)>„(x/3)(/)™(x'7)*(/)fc(x'(5) 

giving the retarded response of the formally non-interacting Kohn-Sham system. The 0^ notation is introduced 
to stress that we deal with retarded quantities. In the above equation (^j(xQ;)'s and e^'s denote respectively KS 
eigenfunctions and corresponding eigenenergies. fj = fri^j), where /T(e) is the Fermi-Dirac distribution function. 
The induced charge and magnetization densities described by the Kohn-Sham susceptibility modify the Hartree and 
exchange-correlation potential, giving rise to a self-consistent problem: the induced densities contribute to the effective 
fields and are, simultaneously, induced by it. The self-consistency is reflected in the second step of the formalism 

X^^(x,x',c.)=xjls(x,x',^)+ II rfxidx2xk's(x,xi,..)te(xi,X2,c.)-l-^^^)x'^(x2,x',..). (9) 



Xi - X2 



The last equation is referred to as "susceptibility Dyson equation", because of its characteristic form, x is the 
density-density response function of the system and describes charge neutral excitations. It is often termed enhanced 
susceptibility. The exchange-correlation kernel, K^c, is defined as a functional derivative of exchange-correlation 
potential with respect to the density 

i^i^[(a(x))](x,x',t-t') = (10) 

evaluated at the ground state values of electronic and magnetic densities. 

The structure of Eq. Q resembles the expression for the susceptibility obtained within the random phase approx- 
imation (RPA) to the many-electron Hamiltoniaui^i^^ However, it is important to note that equation © is exact, 
providing the exact ii'xc were known, whereas RPA corresponds to a specific approximation of the proper polarization 
function. 

The problem of constructing exact exchange-correlation kernel is equivalent to the exact solution of the many-body 
problem and one faces the necessity of approximating it. The most common choice, adopted also in this paper, is 
the ALSDA, in which K^^c is approximated by frequency independent functional derivative of the LSDA exchange- 
correlation potential: 

ifi^[(a(x))](x, X', t-t')^ HS^M^^(, _ _ (11) 

The adiabaticity in this context pertains to the response of the electron gas, which is assumed to be given by the 
instantaneous values of the densities. Furthermore, adopting LSDA implies that the kernel is determined by the 
local value of densities. Recently, there has been a progress in constructing non-local magnetic exchange-correlation 
functionals^^'^'', but their implementation in practical band structure calculations has not been achieved yet. 

From the computational point of view, the main challenge is the evaluation of the unenhanced susceptibility. The 
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direct use of Lehmann representation given by eq. requires the summation over high energy unoccupied states and 
it is practicaUy never used in concrete implementations^, especially for metals. The use of KKR Green's functions in 
the construction of KS susceptibilities avoids this problem. There are several further advantages: it is not necessary 
to consider finite basis corrections and the description of systems without full three-dimensional discrete translational 
symmetry can be rather easily achieved. These issues are discussed deeper in Sec. IIIII 

Not only the calculation but also the analysis of the dynamic susceptibility of a complex system might become 
cumbersome. Because of large number of degrees of freedom the response is governed by a complex overlap of 
various excitation modes.'**' In order to distinguish them, it is convenient to consider the loss matrix, defined as the 
anti-Hermitian part of the susceptibility 

C [x^^ (x, x', + iO+)] ^ i (x^^' (x, x', + 10+) - x''^ (x', X, + iO+)*) . (12) 

Cx has a clear physical interpretation, as the power absorbed from the harmonic driving potential^i 

Si(xt) — Si(x)cosa;i (13) 

reads 

P = J I rfxrfx'S,(x)/:[x'Mx,x',c.)]S,(x'). (14) 

Linear response theory and the fluctuation-dissipation theorem^^iS tell us that the energy absorption signifies the 
presence of excited states of the unperturbed system with the energy lo. Cx is a Hermitian matrix and features real 
eigenvalues and a set of orthogonal eigenvectors. The eigenvectors ^a(x) of £[x*-'(x,x',w)] represent the shapes of 
natural modes of the system. The associated eigenvalues, £[x*-' (x, x',a;)]^, give the rate of energy absorption from 
the external field with the shape fA(x) oscillating with frequency uj. Formally, the number of eigenvalues is infinite, 
but in reality only limited number of them are large, corresponding to physically relevant excitations. The analysis of 
Cx allows their unambiguous identification. The contact with experiment is usually made by evaluating the Fourier 
transformation of the susceptibility as the lii\x^-y^{ci, lo) is probed in the scattering experiments^SiM, cf. App. Rl 



C. Transverse magnetic susceptibility 

Within the non-relativistic LSDA, the Kohn-Sham states of the collinear magnets and paramagnetic systems can be 
characterized by a certain value of the spin-projection. We adopt the convention that the ground-state magnetization 
m(x) points everywhere along the z direction that is selected as the axis of spin-quantization. 

In this case the susceptibility Xks '-'^^^ ^'^^^ independent elements and the following structure 



XKS = 



Xks 


xy 

Xks 








xy 

Xks 


Xks 














Xks 










Oz 

Xks 





(15) 



The response to the transverse (ie. with the direction lying in the xy-plane) magnetic field is transverse and does not 
involve charge density response as opposed to longitudinal magnetic field (along z direction). 

In ALSDA, the exchange-correlation kernel is given by the functional derivative of the LSDA exchange-correlation 
magnetic field 

m 

B,e(m) = B,c(m)— (16) 
m 

with respect to magnetizatio n^^i^^ 

-J = ^^c-R V = [Oij 2~j+~? 2~- vl'J 

ovij orrij m omj m m \ / dm to^ 

The first term gives the response in the direction perpendicular to m = z (the transverse response), while the second 
along the direction of ground state magnetization. We see that the induced transverse magnetization gives rise to an 
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additional effective exchange-correlation magnetic field which is also transverse. Thus, in collinear magnets ALSDA 
decouples magnons and non-spin-flip excitations. 

This is a useful property since it allows us to separate the Dyson equation for the transverse magnetic susceptibility 
from the one for the longitudinal and the charge response. In order to study spin-flip excitations in collinear systems 
it is sufficient to consider transverse components {xx, and xy) of the magnetic susceptibility. It is convenient to 
introduce circular coordinates for magnetization response, magnetic field, and the transverse susceptibility 

m±=m,±im„ B^ = B^±iBy, x*=x""Tix"^. (18) 
They are connected through the following relation 

m=t(x,w) = yrfx'x=^(x,x',w)B±(x',w). (19) 
The explicit form of the KS susceptibility in the collinear system reads 

fcm ^ + ^k - + I" 

Left arrows corresponds to Xks while the right arrows to x^s- 
The Dyson equation for enhanced susceptibility takes the form 

X=^(x,x',w) = Xks(X'^'''^) + y dxiXKs(x,Xi,a;)iv:xc(xi)x=^(xi,x',a;). (21) 

The use of circular coordinates allows us to work with a complex scalar equation instead of a matrix equation in 
Cartesian spin projections. The real exchange-correlation kernel amounts to 

Axc(x) = -MB^^^, (22) 
m(x) 

where _Bxc(x) and m(x) are local values of the exchange-correlation magnetic fields and magnetization density respec- 
tively. They can easily be found once a LSDA calculation of the ground state is performed. 

If one casts the spatial dependence of x^, and Ky^^ in an orthonormal basis, the Dyson equation takes a matrix 
form with the formal solution 

X±(c.) = (I - x^s(^)^-(^))"'xis(^)- (23) 

According to this equation the singularities of the enhanced (physical) susceptibility stem from two sources. The first 
are the singularities of XksC"^) corresponding to electronic transitions between occupied and empty KS states with 
opposite spin projections, so called Stoner transitions: they form the Stoner continuum. The second correspond to 
zeros of the I — x^s^xc- The step from unenhanced susceptibility to the enhanced susceptibility results in a remarkable 
property that the energy absorption can now take place for frequencies outside the Stoner continuum. This signifies 
the new type of the excitations different from the Stoner transitions described by the unenhanced susceptibility. These 
are the spin-waves. Their formation in complex solids and interactions with Stoner excitations are the subject of the 
next section. 



D. Collective modes 



For the sake of subsequent analysis it is convenient to rewrite equation (|23p as 

X"'M-XksH-^-(^). (24) 

We assume that the inverses of KS and enhanced susceptibility matrices exist for the frequencies of interest. The 
following analysis does not require ii'xc to be adiabatic or local. If the Hamiltonian under consideration admits certain 
symmetries, the Hilbert space in which we represent the susceptibility can be decomposed into subspaces where the 
analysis can be performed independent of each other. For example, for periodic solids the quasi-momentum q e f^ez, 
where fisz stands for the first Brillouin zone, is a good quantum number and the susceptibility becomes block diagonal 
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when represented in the Bloch basis. To simphfy notation, q is omitted in this section. 



1. Frequency ljo outside the Stoner continuum 

Away from the Stoner continuum and the singularities of the kernel the matrix x^^(a;o), ujq G M, is Herniitian. 
Let {vx, \mx)~\ be an eigensystem of x^^(wo). A labels eigenvalues. If all vx 7^ 0, the inverse x{^n) is non-singular 
Hermitian matrix and there are no magnetically excited states at loq. However, if there can form a resonance between 
an external magnetic field \B\) and the induced effective exchange-correlation field, such that the conditions 

I^a) =XksK)I™a), (25a) 
\Bx) = KU^o)\mx) (25b) 

are fulfilled, the eigenvalue v\ vanishes. The enhanced susceptibility becomes singular which signifies the formation of 
a collective magnetic excited state (spin- wave) with energy ujq. It is an exact eigenstate of the many-body Hamiltonian 
and its life-time is infinite. 

Let us construct the loss matrix associated with + iO"*") close to ojq in this case. It will allow us to find spatial 
shape of the resonant field and the fluctuating magnetization density associated with the magnon. We expand x^^{z) 
around Wq 

X-\z) = x^'(^o) + {z- luo)Sx-\loo) + 0{{z - c^of). (26) 

We assume that loq is an isolated singular point, i.e. close to wo the matrix x~^iz) is invertible, providing z ^ wo- 
Let Af be the set of indices A corresponding to vx = 0. If there exists more than one vanishing eigenvalue (the case 
of degenerated magnon states), we shall work with linear combinations of corresponding eigenvectors fulfilling the 
condition 

{mx\Sx-\iOa)\my) ^sxSxx', X,X'eM. (27) 
li Sx ^ we can construct the —1st term of the Laurent series of x(z) around ujq 

x(.)=E-^^^+0(l), (28) 

which is sufficient to determine the associated loss matrix 

Cx{oj) = -TrS{uj-ujo)y2 —\mx}{mx\. (29) 

The singularity of x('^o) of the form discussed above signifies that there exists a spin-wave of energy wo which is an 
exact excited eigenstate of the system. Because of the linear character of the response an external field \B) oscillating 
with the frequency cjq will excite the spin-wave and result in a strong energy absorption providing that the spatial 
form of the field and the spatial form of the loss-matrix eigenstate are not orthogonal 

{B\mx) 7^ 0, XeAf. (30) 

The magnetization response to such a field is infinite in the linear response approximation. Physically, condition pop 
determine the fields able to excite mode A. 

A similar argument can be used to prove that the LRDFT features Goldstone mode, i.e. an excitation of vanishing 
energy22iS^ corresponding to a singularity of static enhanced susceptibility. We want to prove that there exists a 
small static transverse magnetic field i?o(x), pointing everywhere in the same direction, say along x-axis, for which 
the magnetization response is infinite. This corresponds to a divergence of susceptibility and the presence of Goldstone 
mode. It is easy to prove thaliS^ 

|5xc)-XKs(0)kGs), (31a) 
|S,e) =i^xc(0)|mGs), (31b) 

where the |Bxc) is a field of the shape of the ground state exchange correlation field, but pointing in the transverse 
direction and, correspondingly, jmcs) has the shape of ground-state value of the magnetization. All the response 
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functions above are static and the expression is valid for general exchange-correlation kernel; it is also fulfilled in 
ALSDA. This is turn means that the response to the transverse static field of the shape of ground-state exchange- 
correlation field is infinite. This proves that LRTDDFT formalism features Goldstone mode. When dealing with 
periodic (uniform) systems, Goldstone mode is necessarily associated with q = 0, because i?xc(x) is a periodic 
(constant) function of x. 



2. Frequency inside the Stoner continuum 



Now we turn to the case that there are Stoner excitations with energy loq. In this case, the resonance condition ((25)) 
cannot be fulfilled exactly anymore, as xks('^o) ceases to be Hermitian and a non-zero phase is introduced between 
the driving field and the Kohn-Sham response. Some eigenvalues of I — Xks('^)-^xc(w) can nevertheless become small 
and strongly enhance certain eigenvalues of Cxitj). The latter might feature a clear peak with maximum at uq. If the 
density of Stoner excitations around wq is small and weakly depends on cj, the peak will take the form of Lorentzian- 
like resonance. Its non-zero width signifies that the corresponding spin-wave is not an exact eigenstate and features 
a finite life time. Physically, the attenuation is interpreted as a consequence of hybridization of the spin-waves with 
Stoner excitations. It is called Landau damping. In the region of dense Stoner continuum, i.e. where the Hermitian 
part of Xks('^) becomes comparable with its anti- Hermitian part, no well defined spin- waves form, leading to the 
phenomenon of spin-wave disappearance. We will visualize these different regimes on concrete examples in Sec. IIVI 

It is important to note that for a particular frequency w both the real part, which determines the magnon energy in 
the first line, and the imaginary part of the KS susceptibility (responsible for damping) are determined by the Stoner 
transitions in the system. The imaginary part is given by "actual transtions" in the sense that only Stoner pairs with 
energy differences equal to ui contribute to it. On the other hand the real part is determined by both the "actual 
transtions" and "virtual transtions", as it ivolves energy integral over all Stoner excitations. 

We remark that the Landau mechanism is the only one leading to the finite life-time of collective excitations when 
adiabatic approximation of the exchange-correlation kernel is invoked. The neglected complex singularities of K^d^jj) 
matrix could lead to the broadening of the spin-wave peaks also outside the Stoner continuum and account for the 
effects like magnon-magnon or magnon-electron scattering. 



E. Sum rules 



The following sum rule®^ provides further insight in the relation between Stoner continuum and formation of spin- 
waves. One shows that the integration of the the loss matrix associated with transverse Kohn-Sham susceptibility 
over all frequencies is related to the ground state magnetization 



dx'C 



= -27rrnGs(x). 



(32) 



The sum rule for x~ differs only by sign in the right-hand side of the above equation. Applying Cauchy integral one 
proves that integrated loss matrix of the enhanced transverse susceptibility is given by exactly the same expression 
providing the exchange-correlation kernel is taken to be frequency independent 



= =F27rmGs(x). 



(33) 



The same relation holds true for the lattice Fourier transformed susceptibilities defined by Eq. (jA5 



dr'C 



iO^ 



= / dr'C 



rfwx=^(r,r',q,( 



=F27rmGs(r) 



(34) 



where r, r' g f^ws- 

The above relation has a clear physical interpretation. The self-consistency condition expressed by the Dyson 
equation leads to the redistribution of the spectral density of spin flip excitations. As long as the frequency dependence 
of the exchange-correlation kernel can be neglected, the integrated spectral weight does not change. However, the 
character of the excitations is now very different including spin-waves and hybrids between spin-waves and Stoner 
transitions. In many cases, in particular for small spin-wave momenta, most of the spectral power is concentrated 
in the spin-wave peaks. The surprising result that the Stoner excitations present in the underlying system cannot 
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be excited by the field of tlie corresponding frequency may be understood as follows. In the energy region of Stoner 
continuum there exists a significant phase shift between the external driving field and the magnetization response given 
by xks due to the large anti-Hermitian part of the latter matrix. This results in the induced exchange-correlation field 
that is out-of-phase with the external driving field. In the next step, the induced exchange-correlation field adds up 
to the external field and produces again a contribution that is out of phase to both external field and initial exchange- 
correlation field. This process to be performed up to the self-consistency leads to the compensation of the external 
field by the contributions of the induced exchange correlation fields and the suppression of the Stoner transitions. 

III. COMPUTER IMPLEMENTATION 

In this Section we shall discuss major issues concerning the actual numerical determination of the response functions. 
The problem can be split into two parts, reflecting the structure of LRTDDFT. First, one finds xks on the basis of the 
knowledge of the KS eigensystem. Second, the Dyson equation is solved in order to determine the enhanced (physical) 
susceptibility. 

The determination of the KS susceptibility is the most computationally demanding part of the calculations. We 
begin with formally non-interacting KS Green's function 

G.,(x,x-,z)^^ ^^(""^^^("^^)\ zeC, (35) 

Z — Ej 
j ■' 

where (/)j(xa)'s and e^-'s denote KS eigenfunctions and corresponding eigenenergies. The numerical evaluation of G is 
based on Korringa-Kohn-Rostoker method. The actual representation of the GF in the multiple scattering formalism 
is given in App. [Cl 

By applying Cauchy theorem one obtains the following expression for the KS susceptibility^^'^^ (7 > 7' = 0^)- 

1 /"^ 

X^g(x,x',a; + i7) = -— y de x 

(/T(e + i7')5y (x, x', e + + i7 + i7', e + (7') 
-/T(e - i7')S'ij (x, x', e + w + i7 - (7', e - [7') 
+/T(e + i7')'5»j (x, x', e + 17', e - - [7 + [7') 

-/T(e - iy)^.^ (x, x', e - i7', e - - [7 _ [7')) (36) 
where the S is the product of two KKRGFs 

S'.y(x,x',Zi,Z2) EE Cr^^CT:^^G'/3^(x,x',Zi)G5a(x',X, Z2)- (37) 

When working in k-space, the product above is transformed into a convolution of two Green's functions over the 
Brillouin zone, cf. App. [Cl The convolution converges badly for energies zi.2 close to the real KS energies. This 
in turn precludes the direct determination of the retarded KS susceptibility."^^ For computational expediency one 
calculates the susceptibility for a set of the points in the upper complex semi-plane, away from the real axis. As 
shown in App. [B] in this case, the necessary energy integrations in eq. (|36p can be reduced to integrations along finite 
complex contours, away from the poles of the GFs. At the end of the calculations the resultant susceptibility can be 
analytically continued to the real axis to recover the real time dynamics. 

Note the presence of the Fermi-Dirac weight in eq. ([36]) . The formulation presented avoids integration over energies 
above Fermi level. It is so, because KKRGF contains information about all KS states, cf. eq. Thus, the serious 

technical problem of including explicitly unoccupied KS orbitals, as it is the case when evaluating xks directly from 
eq. ([8|), is avoided.^ No finite basis set corrections are necessary, either, when working with KKRGF— i22r^, as 
the multiple scattering problem is solved separately for every complex energy. In this context, it is not clear how 
large is the error introduced by the minimal energy independent d-symmetry basis used by Lounis et aL ^^i^^ , an 
approximation avoided in this work. In our current study the major advantage of using KKRGF is the possibility of 
an efficient description of systems featuring surfaces and interfaces^^^, in particular films absorbed on a substrate. 
The construction of xks taking into account specific representation of KKRGF is given in App. [Cj 

It is convenient to cast the (x, x') dependence into a separable basis when solving the Dyson equation. This aspect 
is discussed in Appendix [Xj Subsequently, the equation (j2ip can be solved by matrix inversion or, for frequencies 
away from the spin-wave poles, iteratively. The CPU time necessary to solve the susceptibility Dyson equation is 
negligible comparing to the time of computing the KS susceptibility. 
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Below we turn to the important question of the consistent description of the Goldstone mode. As discussed in Sec. 
IIIDl x(q = 0, w = 0) is a singular matrix with a diverging eigenvalue signifying the formation of the Goldstone mode. 
Seen alternatively, cf. eqs. pip , the matrix 

I?=l-xis(0)X.c (38) 

features a vanishing eigenvalue corresponding to the eigenvector representing ground state magnetization |togs)- This 
is the consequence of the spin rotational invariance of the problem. 

Because of the numerical inaccuracies the condition of zero eigenvalue in the limit of q — is not satisfied exactly 
since the calculated Kohn-Sham susceptibility does not correspond exactly to the ground-state exchange-correlation 
magnetic field and ground state magnetization. Below we discuss the method to compensate for this error. Numerical 
diagonalization of matrix 2? gives one eigenvalue that is very close to zero and much smaller than other eigenvalues. 
In our calculations for systems discussed in this paper this eigenvalue is typically of the order of 10^'^ whereas other 
eigenvalues are of the order of unity. It is easy to verify that this small eigenvalue corresponds to the eigenvector 
very close to |mGs)- To obtain the matrix with zero eigenvalue we proceed as follows. Upon the diagonalization of 
the calculated I?, the small eigenvalue is set to zero and the remaining eigenvalues are left unchanged. Using the 
original eigenvectors of T) the corrected diagonal matrix is transformed back into non-diagonal form. Because of the 
high precision of the diagonalization routines the corrected matrix I?corr will have a zero eigenvalue with a very high 
accuracy. I^corr is now used to find the corrected exchange-correlation kernel 

^xr-(xis(0))"'(l-2') (39) 

This corrected kernel corresponds to the calculated Kohn-Sham susceptibility Xks(^) ™^ ^^'^ sense of the fulfillment of 
the Goldstone theorem and is used for the calculation of the enhanced susceptibility for all q vectors and frequencies. 
For the example of bcc Fe the uncorrected small eigenvalue of T> matrix in our calculation was e w 3.6 • 10"'^ which 
corresponded to the shift of the Goldstone mode from zero by about 2.3 meV (compare to the energy width of magnon 
band of order of 0.5 eV). 

The methods to correct for the numerical deviations from the requirement of Goldstone theorem suggested by 
Lounis et al^, Kotani et al^ and §a§ioglu et al^^ are all based on the same general idea of bringing in necessary 
correspondence of the unenhanced susceptibility and underlying potential expressed by equation pif) . That is also 
the essence of our approach. Our correction of one eigenvalue of order of 10"'^ should be compared with 2 — 12% 
correction for the exchange-correlation kernel reported by Lounis et al.^^ and 50% correction for the screend Coulomb 
interaction of §a§ioglu et alJ^ that demonstrates the robustness and accuracy of our method. 

Finally, we mention that the computations of KS susceptibility can be massively parallelized and we resorted to the 
Message Passing Interface (MPI)^'*"^^. The calculations of non-enhanced susceptibility for every frequency and wave 
vector are essentially independent from each other and can be performed on separate processors. No inter-processor 
communication is involved in this mode and the effectiveness of parallelization increases linearly with the number of 
processors even for their number as big as several thousands. The most time consuming calculations presented in this 
paper required several days of computational time on a modern 64 processor machine. 

IV. NUMERICAL RESULTS AND DISCUSSION 

The standard picture of magnons in itinerant ferromagnets is based on the random phase approximation treatment 
of the uniform electron gas.""^ In this model the only spin-wave branch evolves from the zero energy Goldstone mode. 
At g =: 0, the whole spectral power of the Stoner continuum is transferred to the magnon pole, cf. Sec. Ill Dl As 
the momentum increases the energy of the magnon rises (~ for small momenta) and its life-time remains infinite 
until the magnon band makes a contact with the Stoner continuum, where the spin-wave abruptly looses its resonant 
character and cannot be considered as a well defined excitation. The strength of the latter effect, called "spin- wave 
disappearance"—, follows from the properties of the uniform electron gas. For the frequencies corresponding to Stoner 
transitions the imaginary part of the KS susceptibility is practically everywhere comparable in magnitude to its real 
part. 

Real materials feature much richer spin dynamics. There are at least two important reasons for this qualitative 
difference with the uniform electron gas. The first is the formation of atomic magnetic moments resulting from highly 
non-uniform spatial distribution of charge and spin densities. Although the atomic magnetic moments are formed 
by the itinerant electrons, the strong intraatomic exchange interaction keeps them well defined even at elevated 
temperatures. This feature is the reason of the usefulness of the Heisenberg model of interacting atomic spins in the 
discussion of the physics of itinerant-electron magnets. The strongly non-uniform distribution of magnetization on 
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the atomic length-scale leads to two further effects not present in the uniform electron gas. First, multiple spin-wave 
branches form for momenta from the first Brillouin zone of the systemi--'^- Second, not the entire spectral weight of 
the Stoner continuum is shifted to the Goldstone mode for q — 0; part of it resides in the optical spin-wave modes 
and the residual Stoner continuum. 

The second reason is the complexity of Stoner continuum reflecting the rich electronic structure of real materials. 
The latter is always of a multiple-band type; the bands differ strongly from each other in the character of the 
hybridization of the atomic s, p and d states of the same and different atoms. It is common to characterize the 
electronic structure of a ferromagnetic system by the exchange splitting A = Im where / is the so-called Stoner 
parameter of the material and m is the magnetic moment per atom. It is, however, very important that A characterizes 
the energy splitting only between the Sd states of similar spatial form. On the other hand, the Stoner continuum 
contains the spin-fiip transitions between all available pairs of electron states with opposite spin projections. For a 
given energy and wave vector the value of the spectral density of the Stoner transitions is determined by the number 
of the states available for the transitions weighted with the respective matrix element. It depends crucially on the 
overlap of the wave functions of the initial and final electronic states. For the states separated by the exchange 
splitting A, the overlap of the wave functions is large which leads to high spectral density of the Stoner transitions 
around this energy. As follows from the discussion in Sec. IIIDl a strong destructive influence of the Stoner transitions 
on the spin-wave states is expected for spin-waves forming in this frequency range. In many systems, however, the 
value of A exceeds substantially the spin-wave energies and the corresponding Stoner transitions cannot contribute 
to the Landau damping. The latter is governed by the low-energy part of the Stoner spectrum involving the electron 
transitions between states with possibly very different orbital character and therefore with a low transition probability. 
The low intensity of the Stoner spectral density at given q and lu leads to the weak Landau damping. As a result, 
the spin-waves corresponding to these q and uj remain well-defined. 

Obviously, the form of the Stoner continuum, and in particular its low-energy part, is strongly system dependent 
reflecting details of the electronic hybridizations in the system. Not only the chemical composition but also the 
dimensionality of the material exercises strong influence on these properties. In ultrathin films, the presence of 
the nonmagnetic substrate brings the states of the magnetic film in contact with the states of the substrate, which 
complicates further electron spectrum and correspondingly Stoner continuum. In Ref.*"'^ we showed that arguments 
based on the simplified analysis of the dimensionality aspects are not sufficient to predict the properties of the spin- 
wave damping in the magnetic films. The concrete details of the electronic structure are essential. Depending on 
these details, the interpretation of the features of the spectral density of spin-fiip excitations varies from the infinitely 
living Heisenberg-type spin-wave modes, through moderately damped well defined magnons up to the spin-wave 
disappearance effect. 

Below we report a systematic comparison of the spin-wave damping in the bulk and film forms of Fe, Co and Ni. We 
show that the trends in the dimensionality dependence of the spin- wave properties are different for different elements 
and interpret them on the basis of the analysis of the underlying electronic structure and spin-flip spectrum. We 
primarily analyze spin- wave spectra based on the eigensystem of the loss matrices. The energy uig of the magnon is 
identifled with the energy position of the spin-wave peak maximum. The full width at the half maximum (FWHM) 
of the peak gives the inverse life-time of the excitation. We present also the cases where the energy dependence of the 
eigenvalues cannot be described by a well deflned peak. We briefly assess the applicability of the adiabatic methods 
to estimate frequencies of short wave-length magnons. 

A. Fe 

1. bulk bcc Fe 

According to our calculations, the spin-wave disappearance effect is particularly pronounced in the bcc Fe, cf. 
Figs. [T] and ^ Above the critical energy of 82 meV corresponding to qc = 0.35^ wave-vector the spin-wave band 
enters abruptly a region of dense Stoner continuum and the intensity contained in the spin-wave peak drops by an 
order of magnitude. The strong loss of spin- wave intensity above around 80 meV was clearly observed in the neutron 
scattering experiments, cf. e.g. Refi^. The effect is rather anisotropic and sets in mainly along F-H direction, which 
nicely matches the experiment as welli^ Fig. [T] shows clearly that the strong damping is directly related to the 
appearance of low energy Stoner excitations at qc- The origin the these low energy Stoner transitions can be traced 
back to relatively large density of majority d spin electrons at the Fermi level in this material. 

Also the energies of spin-waves in the adiabatic local density approximation correlate nicely with experimental 
values. ^'^^^ According to our data for q < qc the dispersion can be well represented by the biquadratic flt 



UJo{q)^Dq^{l-jq^). 



(40) 
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FIG. 1: Examples of spin-flip spectra in bcc iron, atomic units, for different momenta along (100) direction. The largest 
eigenvalues of respectively enhanced and KS loss matrix are shown. Low energy spin-wave peaks have simple Lorentzian shapes 
small line- widths and carry substantial spectral power. Above the critical energy of 82 meV corresponding to qc = 0.352^ 
wave-vector the spin-wave band enters abruptly a region of dense Stoner continuum. The spin-wave peaks become broad, 
acquire irregular non-Lorentzian shapes reflecting the energy dependence of the density of Stoner continuum, and carry much 
smaller spectral weight. 



with parameters D = 252 meVA and 7 = 0.28 A . Experimentally reported values for D vary between 266 meVA 
and 307 meVA'^>27-ioo 

For the well defined, low energy spin-waves, the loss matrix of enhanced susceptibility features only one large 
eigenvalue. The associated eigenvector corresponds to the practically rigid rotation of magnetic moments. This 
results justifies the applicability of Heisenberg model to describe low energy magnons. For high energy strongly 
damped magnons, the spectrum contains only one dominating eigenvalue as well, but the corresponding eigenvectors 
involves small, but clear deformation of the atomic moment. On the other hand, the unenhanced susceptibility features 
multiple eigenvalues of comparable magnitude, cf. Fig. [31 corresponding to transitions between different KS orbitals. 
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FIG. 2: Spin-waves of bcc iron obtained using LRTDDFT implementation described in this paper. Solid circles (•) correspond 
to the maximum of spin-wave peak, while the error bars denote full width at half maximum (FWHM) of the peak. Solid line 
denotes spin- wave energies obtained from magnetic force theorem (MFT). Strong increase of Landau damping in all directions 
in the Brillouin zone (spin- wave disappearance effect) is seen for spin- wave energies above 82 meV. High energy spectrum along 
FHN directions is dominated by Stoner excitations and cannot be represented by a simple Lorentzian. Representative examples 
of the spectral power functions are given in Figure [T] 




Lj [eV] 



FIG. 3: The largest eigenvalues of loss matrix associated with unenhanced (KS) susceptibility for a selected momentum in bulk 
bcc Fe. There are several eigenvalues of comparable magnitude corresponding to transitions between different KS orbitals. On 
the contrary, the loss matrix associated with the enhanced susceptibility features only one eigenvalue, identified as a single 
spin-wave mode of the monatomic material. 



2. Fe films 

Contrary to the bcc bulk case, the free standing monolayers of Fe, both of (100) and (110) crystallographic orienta- 
tions, feature well defined spin- waves for all momentap^i^"— see also Figs. 3] and [S] This property can be traced 
back to the enhanced exchange splitting of d-symmetry orbitals and the band narrowing in the free film. It effectively 
removes the majority-spin d states from the Fermi level— resulting in the small density of low-energy Stoner 
excitations. The same behavior is observed in unsupported Ni and Co (100) monolayers, which will be discussed later. 

When the monolayers of Fe are deposited on a substrates the spin-dynamics is modified. We demonstrated recently^ 
that the impact of the substrate varies strongly depending on the substrate and its orientation. Explicit LRTDDFT 
calculations for Fe/Cu(100) system have been performed by us for one**^ and three monolayers'^ coverages. The 
damping increases somewhat compared to the free standing case, but the spin-waves are well defined in the whole 
two-dimensional Brillouin zone, cf. Fig. 21 The substrate-induced renormalization of magnon energies turns out to be 
small. In the case of Fe/Cu(100) the states of the substrate hybridize rather weakly with the states of the magnetic 
overlayer. The films feature small density of majority-spin states at the Fermi level and the spin dynamics of the films 
resembles qualitatively the one of free standing layers. 

The situation is dramatically different in Fe/W(110), as evident from Fig. [SJ The substrate renormalizes magnon 
energies and strongly enhances the damping. Here, an important role is played by the interface electronic complexes, 
formed by hybridized surface states of W(llO) and electron states of the film^. The complexes provide an efficient 
source of Stoner pairs in the region of magnetic overlayer. The strongly damped collective precession has been 



15 




FIG. 4: Spin-waves of Fe(lOO) monolayer free-standing and supported on Cu(lOO) surface. Solid lines correspond to cJo(q), 
while the width of the shaded region to FWHM. The spin dynamics of the free standing and supported monolayers differ weakly. 
From^. 




FIG. 5: Spin-waves of Fe(llO) monolayer free-standing and supported on W(llO) surface. Solid lines correspond to ajo(q), 
while the width of the shaded region to FWHM. The surface state of W(llO) leads to a qualitative change in the spin dynamics 
of the absorbed film. From*^ . 



observed experimentally at the zone boundary in 1 ML Fe/W(110) systeni^°^. We emphasize, however, that even for 
Fe/W(110) the region in the Brillouin zone featuring the spin-wave disappearance effect is small compared to the 
bulk case. Experiments^ shows also that the damping-to-energy ratio in 2 ML Fe/W(110) film decreases compared 
to the single monolayer case. The trend correlates well with our calculations, cf. Fig. [HI We remark that the system is 
characterized by quite complex structure and in our calculations we took into account the atomic relaxationsi2^-i2^. 

In 2 ML Fe/W(110) system there are two types of non-equivalent Fe atoms and the loss matrix associated with 
the enhanced susceptibility features two large eigenvalues corresponding to the acoustic and optical spin- wave modes. 
The spectra are qualitatively very similar to those of hep Co presented in Sec. IIVB II 

The observed magnon energies are roughly 40% smaller than predicted by our theory. The reason of the discrepancy 
is still not clear, especially in the light of good performance of ALSDA in the bulk bcc iror>i2S. It is worth to remark 
here that the LRTDDFT performs much better compared to model Hamiltonians, in which the spin- wave energies 
of 1 ML Fe/W(110) are grossly overestimated^^i^iiiiSi. Furthermore, the experiment did not reveal any presence of 
an optical branch in the 2 ML case. It has been conjectured^^" that SPEELS could probe only the modes with 
significant amplitude in the top layer, because of a limited penetration depth of electrons. However, our calculations 
do not predict the formation of modes localized at the surface of the film. Therefore, the excitation of both types of 
spin-waves should be expected. 

For large wave vectors the both spin-wave peaks in the spectral density are substantially broadened and might 
appear as a single feature in the spectrum, especially if the finite SPEELS resolution is taken into account. On 
the other hand, in the center of the zone the optical mode should be discernible. Possibly, the optical mode, being 
substantially broadened even for g = 0, is lost in the background of the signal dominated by the acoustic mode. 

We remark that in the case of 1 ML Fe/W(110) the MFT yields the spin-wave energies*!^ii2^ of similar values as 
LRTDDFT; the same is the case for 1 ML Fe/Cu(100) studied by Pajda et al^. Udvardi et al^ considered also 
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FIG. 6: Spin-waves of 2 ML Fe/W(110). Two magnon modes can be discerned in the whole two dimensional Brillouin zone. 
Experimental points come from Ref.iiS^. 



the relativistic corrections. In the relativistic limit the spin is coupled to the lattice and the Hamiltonian looses its 
spin-rotational invariance. As a consequence the lowest energy magnon does not have vanishing frequency and the 
Goldstone mode disappears. Udvardi et al. show that in light transition metals the relativistic corrections to the spin- 
wave spectra are of minor importance on the energy scale of exchange interactions. Simultaneously, they predict the 
appearance of a series of weakly dispersive spin- wave bands associated with the dynamics of small magnetic moments 
induced in the Cu substrate. These additional resonances are absent in our data, also in the case of Fe/W(110), 
where the moment induced in the interface layer of W is sizable (ss 0.14 /is) and antiferromagnetically aligned with 
the film magnetization. As it has been already discussed^^^ the appearance of additional spin-wave modes related to 
the induced moments is an artifact of adiabatic approach where small induced moments are treated as independent 
adiabatic degrees of freedom. The methods based on the evaluation of magnetic susceptibility are better suited for 
the investigation of such "induced" magnetization dynamics^S.. 



B. Co 



At ambient pressure and low temperature cobalt is a ferromagnet featuring e (hep) structure. -i*!^ As the temperature 
increases a transition to 7 (fee) structure occurs at around 750 K and above Curie temperature of around 1400 K the 
systems becomes paramagnetic. To our knowledge, all previous ah initio studies of spin-waves in Co employ adiabatic 
approximation i^'^^1^^ 
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c) ^[cV] d) '-'[cV] 

FIG. 7: Examples of spin-flip excitation spectrum in hep Co, atomic units, for two momenta along FA direction. Panels (a) and 
(b) present the two dominant eigenvalues of the loss matrix, Cx{'ii^)\ the eigenvalues correspond to the acoustic and optical 
spin- wave branches. As expected from symmetry arguments, they become degenerate for Qz — ^. Panels (c) and (d) show the 
imaginary part of the Fourier transformed susceptibility, lmxKii{ci, (jj) , for momenta inside and outside the first Brillouin zone. 
The momenta differ by a reciprocal lattice vector. ImxKK(q, (jj) is probed in the neutron scattering experiments. By varying 
the momentum one can access the optical and acoustic spin wave branch. 



1. e(hcp)-Co 

Experimental lattice constant of hep cobalt Ohcp = 4.738 au was used.^^^ The calculated magnetic moment per Co 
atom reads 1.61 /zb- 

The loss matrix associated with the susceptibility features two large eigenvalues for q in the first Brillouin zone and 
we interpret them as acoustic and optical spin-wave branches, cf. Fig. [T] The positions and widths of the peaks are 
presented in Fig. [8] together with adiabatic spectra obtained from magnetic force theorem. Our method, MFT and 
STSM yield rather similar energies. For = ^ plane, the both modes are degenerate; it is the consequence of the 
symmetry of the hep lattice. Theoretical spin-wave energies correspond very well to the experimental results along 
FM direction, but are larger along c-axis (FA) direction, cf. Fig. [S] For small momentum transfers the peak position 
can be very well described by the biquadratic fit, cf. eq. |40l Table H] presents the results of the fit for q < 0.3^^ 
along different directions. Note that for the hep system the dispersion relation for small q is isotropic in the basal 
plane. Our parameters D weakly depend on the direction in the reciprocal space because of the limited accuracy of 
the biquadratic function in the fitted interval of the wave vectors. 



direction D 


[ meVA'] 


7[A^] 


FM 


539 


0.27 


FK 


520 


0.26 


FA 


529 


0.32 



TABLE I: Parameters of the biquadratic fit of Eq. (I40p for different direction in Brillouin zone in hep Co. 
Additionally to the spectrum of the loss matrix, £[x(q, w)], we present in the Fig.[7]the imaginary part of the Fourier 
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FIG. 8: Spin-waves of of hep cobalt, (a) Position of the peak uoo{<l) obtained from the enhanced susceptibility. The energies of 
spin-waves are degenerated within the numerical error along ALHA directions. The experimental energies come from Refs.— li^^. 
(b) Adiabatic magnon spectra obtained by means of magnetic force theorem, (c) half-widths at half maximum, FWHM= 2/3(q). 



transform of the enhanced susceptibility, ImxGG(q, for two selected momenta. By varying the momentum transfer 
in the inelastic scattering experiments, one can modify the intensity of peaks coming from the optical and acoustic 
spin-waves. ^''^^ In the example of hep Co, along the FA direction, only the acoustic magnons are observed in the first 
Brillouin zone, while the optical ones are detectable only in the second Brillouin zone. Perring et alr^ succeeded to 
detect spin-waves for q beyond the first Brillouin zone along (00^) direction, thus accessing the optical SW branch 
of FA segment. The formalism of this paper predicts that the two branches exist in the whole Brillouin zone. The 
acoustic mode corresponds to the moments oscillating in phase and the optical mode in anti-phase. The energy of 
the optical mode is overestimated in ALSDA. The optical magnons are of much shorter life-time than the acoustic 
magnons because of the higher density of the Stoner transitions, as expected for higher energies. Characteristic peaks 
in the FWHM curves correspond to the areas where the spin- wave branch crosses the region of larger density of Stoner 
excitations. The damping of spin-waves is moderate in the case of Co and peaks have well defined Lorentzian shape. 
All the majority spin d electrons are occupied and located rather far from the Fermi level. As a consequence, the 
low energy Stoner excitations involve primarily and to d'^ transitions which results in the small Stoner intensity 
because of weak overlap of the wave-functions of the initial and final states. 

Earlier theoretical works on spin susceptibility of hep cobalt utilize empirical tight binding schemsii^iii^. They 
predict correctly the energies of acoustic modes, yielding however too low values of optical modes, in fact in a range 
where they were not detected despite being experimentally observable. Except for FA segment the optical branch 
energies predicted in our study lie above 0.5 eV, the maximal energy addressed in the calculations mentioned and one 
guesses that certain complex structure of the spectral density of spin-flip excitations (cf. the case of bcc Fe in Fig. [T]) 
might have been erroneously identified as the optical mode. 
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FIG. 9: Spin-waves of fee Co. Solid cireles (•) eorrespond to ajo(q), while the error bars denote full width at half maximum 
(FWHM) of the peak. Solid line denotes spin- wave energies obtained from MFT. Solid triangles (a) stand for the experimental 
estimation of the dispersion relation by Balashov.-^^S 

2. -yffccJ-Co 

The lattice constant employed in the calculations was determined from the condition of equal atomic volumes for 
hep and fee systems, yielding Ofcc = v^ahcp- Ground state magnetic moment is very close to the one obtained for 
hep phase. 

The loss matrix of the enhanced susceptibility features only one large eigenvalue. The high wave- vectors magnons 
show slightly higher energy, when compared to the MFT results, cf. Fig. [S] The spin-waves are clearly damped, but 
the spectrum can be well described by a Lorentzian peak for all momenta; no spin-wave disappearance is observed. 

Inelastic neutron scattering experiments revealed magnons of slightly smaller energies J^ii^S, The stiffness con- 
stant estimated from our data amounts to D = 492 meVA , while D = 369 meVA and D — 356 meVA 
were reported. The results of spin polarized electron energy loss spectroscopy^ and inelastic scanning tunneling 
spectroscopy-^*^ match rather well the neutron scattering data. The difference with our calculations can be traced 
back to the finite temperature of the experiment. For the data extracted from the film measurement o^^i^^" the Cu 
substrate might influence the results, but our calculations reported in the next section exclude this possibility. In 
the case of bulk (neutron) measurements fee Co must be alloyed with about 6% of Fe for the sake of stability, which 
might further contribute to the quantative differences between the theory and experiment. 

Rather limited data exists regarding the life-time of magnons. Vollmer et al^ were able to perform constant q 
scans and provided an estimate of FHWM about 40 — 75 meV. This matches quite well our low energy results. 

3. IML Co/Cu(100) 

The spin-waves of 1 ML Co/Cu(100) are clearly of higher energy than in the bulk fee case, see Fig.[TUl Interestingly, 
the Landau damping is slightly smaller. Compared to the free standing monolayer we observe almost no magnon 
energy renormalization and moderate increase of the Landau damping. This behavior is qualitatively very similar to 
the IML Fe/Cu(100) case. 

C. Ni 

1. bulk fee Ni 

The transverse magnetic susceptibility of bulk fee Ni has been earlier studied using different first-principle 
approache9^2i^"2£. All previous calculations and the present study performed by us, cf. Fig. [11] yield a very similar 
picture. The spin dynamics of nickel differs strongly from that of iron, in particular, no spin- wave dissapearance is 
found. Magnon peaks are well defined in the practically entire Brillouin zone, characterized by a relatively small 
Landau damping. This is rather surprising, realizing that the Stoner continuum is pronounced in Ni at relatively low 
energies, corresponding to the exchange splitting predicted by LSDA to be only about 0.7 eV. The effect of the low 
damping of the Ni magnons has its roots in the fact that the spectral power of the Stoner continuum is pronounced 
mostly around the exchange splitting energy and is not spread over a wide range of energies as in Fe. Furthermore, all 
the majority spin d electrons are occupied and located well below the Fermi level. Therefore, the Stoner transitions 



20 



> 

CD 



3 



1.2 
1.0 
0.8 
0.6 
0.4 
0.2 
0.0 





1A/TT Pri/'inn'lfroo 














r X M r 



> 

CD 



f 



1.2 
1.0 
0.8 
0.6 
0.4 
0.2 
0.0 





^ 1MT, Cn/C^^(^(^C\^ 







X 



M 



FIG. 10: Spin-waves of Co(lOO) monolayer free and supported on Cu(lOO) surface. Solid circles (•) correspond to ajo(q), while 
the error bars denote full width at half maximum (FWHM) of the peak. Solid line denotes spin-wave energies obtained from 
MFT. 
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FIG. 11: Spin-waves of fee nickel. Solid circles (•) correspond to a;o(q), while the error bars denote full width at half maximum 
(FWHM) of the peak. Solid line denotes spin-wave energies obtained using MFT and STSM; the two methods differ significantly 
for this system. Most of the spectrum can be described using Lorentzian distribution function, except most energetic spin-waves 
along [^^0] direction. An example of such a irregular spectrum is presented in Figure [121 

involving these states are not effective at the energies characteristic to magnon excitations in fee Ni. In this respect 
the system is similar to hep and fee Co. 

Another special property of Ni is a pronounced non-monotonous dependence of damping on q, seen most clearly 
along [^00] in the Brillouin zone, cf. Fig. 1131 At the X point the spin- waves are very weakly Landau damped. This 
behavior originates again from the very narrow Stoner continuum at the X point; it is almost entirely concentrated 
around the exchange splitting energy, while at A = [5OO] point it has a significant contribution at smaller energies. 
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FIG. 12: An example of spin-flip excitation spectrum in Ni, — £[x(q, a;)];^, atomic units. Only eigenvalue of the loss matrix 
dominates. One clearly sees the energy dependence of the Stoner continuum reflected in the behavior of enhanced susceptibility. 
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FIG. 13: Two examples of spin-flip exeitation speetrum in Ni, atomic units, along (100) direction. The non- monotonous 
dependence of the damping with the magnon momentum is clearly seen. Additionally, close to the X point in the Brillouin 
zone a clear coexistence of the spin-wave peak and a Stoner continuum peak can be observed. 
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FIG. 14: (a) Double peak feature in bulk fee Ni. (b) Behavior of the corresponding Kohn-Sham susceptibility. 



This leads to a very interesting effect at the X point where the spin-wave peak (identified by a small eigenvalue of 
I — X^c^{0)K:x^c matrix) appears just below a peak of the Stoner continuum. Existing experimental studies point to a 
rather monotonous increase of the magnon liiiewidth with increasing momentaii^^ii^l. However, the momenta close to 
the zone boundary are difficult to be probed experimentally using thermal neutrons. We conjecture that the coexisting 
magnon and Stoner peaks were not separately resolved giving rise to one broad spectral feature. 

Much attention has been paid to the appearance of so called "optical spin-wave mode" experimentally detected in 
fee Ni22iii21. In our calculations a low energy double magnon peak structure appears only in a very small part of 
the Brillouin zone around q = (0.15, 0, 0)27r/a, cf. Fig. [Hk. The feature is strictly absent in the Heisenberg model, 
as it features only single degree of freedom in monatomic Ni. We note, however, that (i) the loss matrix of the 
enhanced susceptibility features only one large eigenvalue in this energy region, (ii) both peaks are associated with 
similar eigenvectors and therefore describing similar magnetization shape. In this respect the double peak structure 
clearly differs from the "real" optical spin-wave mode of hep Co arising from the presence of two magnetic atoms 
in the primitive cell. The splitting of the peak was explained by Karlsson et al^ by analyzing the non-enhanced 
susceptibility in the region of the double peak structure, cf. Fig. 114b . characterized by an appearance of rather narrow 
low energy Stoner peak. The corresponding real part of the susceptibility changes non-monotonously with the energy 
bringing an eigenvalue of I — Xks(^)^^(^ close to zero twice on a the energy scale of the peak splitting. 

Unfortunately, the agreement with experiment is much worse regarding the magnon energies. The stiffness con- 
stant extracted from our data equals 851 meVA and is roughly twice as large as the experimentally observed 
374 — 433 meVA -'-^^i-'-^^ . At larger momenta, the same factor two discrepancy between measured and calculated 
magnon energies is seen. The problem is associated with the overestimation of the exchange splitting in LDA for 
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FIG. 16: Spin-waves of Ni(lOO) monolayer free and supported on Cu(lOO) surface. Solid circles (•) correspond to a;o(q), while 
the error bars denote FWHM. Solid line denotes spin-wave energies obtained from MFT. Spin-waves with large momenta are 
not well defined excitations in the case of the absorbed film, while for the free standing film they exist in the whole Brillouin 
zone. 



The adiabatic STSM and the results based on the dynamic susceptibility calculations, agree very well with each 
other but yield clearly higher energies than MFT, cf. Fig. [TT] As it has already been pointed out by Grotheer et al^, 
the systematic error of the MFT method'^^ is particularly pronounced in Ni, owing to its small exchange splitting. 
Nevertheless, the spin- wave stiffness obtained from all adiabatic methods is identical to the one obtained from dynamic 
susceptibility, cf. FigfT5l 



2. 1 ML m/Cu(100) 



Similar to the case of Fe and Co, the magnons in free standing monolayer of Ni (100) are weakly damped, as seen 
in Fig. 1161 Their energies are much higher than in the fee Ni bulk. However, upon absorption on Cu(lOO) surface, the 
spin-dynamics of the systems changes dramatically and in the way that is in strong contrast to the other two transition 
metals. Spin waves in 1 ML Ni/Cu(100) are defined only close to the center of the Brillouin zone and strongly damped 
for larger momenta. Two examples of spin-flip spectra are shown in Fig. [T71 The hybridization of the states of the Ni 
overlayer and the Cu substrate states leads to (i) reduction of the exchange splitting (the Stoner excitations for g = 
are centered around 218 meV vs. 700 meV in the bulk), reflected in the smaller magnetic moment (0.27 /ib vs. 0.62 iib 
in the bulk) , (ii) larger energy width of the Stoner continuum. This two effects result in the enhanced density of low 
energy Stoner excitations, washing out most of sharp high energy spin-wave features. In IML Ni/Cu(100) localized 
moment picture of spin excitations (i.e. Heisenberg model) fails altogether. The Landau maps of IML Ni/Cu(100), 
cf. Fig. [in] are qualitatively similar to the to case of the Fe and Co monolayer*^ but characterized by much broader 
hot-spots and more intense diffused background. This property reflects strong hybridization of the electronic states 
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FIG. 17: Example ol spectra in 1 ML Ni(lOO) and 1 ML Ni/Cu(100), atomic units, q = (0.5, 0)27r/a corresponds to X point in 
the Brillouin zone. In the supported monolayer, for large wave-vectors, the Stoner continuum becomes very broad (in contrast 
to the free standing film) and the Landau damping completely washes out any sharp spin-wave features. 




FIG. 18: Landau maps of 1 ML Ni(lOO) and 1 ML Ni/Cu(100): intensity of Stoner transitions with momentum q — (0.5, 0)27r/a 
(X point) and energy cjo = 400 meV in the Ni layer resolved for different final k-vectors in the first Brillouin zone. The Stoner 
states cause the damping of magnons presented in Fig. [171 

in the film with the substrate electrons. 



V. CONCLUDING REMARKS 

The spin- flip dynamics of elementary Sd transition-metal ferromagnets is rich and strongly system dependent. The 
magnons of Co live relatively long for momenta in the whole Brillouin zone. One might be tempted to associate this 
property with the large exchange splitting of bands, which usually results in the Stoner continuum pronounced at 
high energies, but such explanation fails badly for Fe, which features even larger exchange splitting and at the same 
time severe Landau damping of spin- waves. The effect is even more spectacular in Ni, where the band splitting is only 
a bit larger than the typical magnon energy. Owing to the compactness of the Stoner spectrum and the small density 
of majority spin d states at the Fermi level, Ni features rather long-living magnons in most of its Brillouin zone. We 
see that these are fine properties of Stoner continuum which determine the spin- wave attenuation. Obviously, the first 
principle approaches based on the calculation of transverse magnetic susceptibility are indispensable in the consistent 
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description of spin dynamics in real materials. 

This statement holds particularly true for the spin excitations of ultrathin films. Free standing monolayers studied 
in this paper shows usually well defined long-living spin-waves for all momenta. Upon absorption of the magnetic 
film on a nonmagnetic substrate the damping generally increases, but the details of the spin-flip dynamics are very 
sensitive to the details of the electronic hybridization between the electron states of substrate and film. Monolayers of 
Fe (supported both on Cu(lOO) and W(llO) surfaces) feature relatively long-living, well defined spin-waves in most of 
the two-dimensional Brillouin zone, in striking contrast to the bulk bcc phase. In 1 ML Ni/Cu(100) spin- waves exist 
only for small momenta, the spin-flip spectrum in the rest of the Brillouin zone is dominated by incoherent Stoner 
excitations. The situation is again opposite to the one in the bulk fee Ni bulk characterized by generally long-living 
high energy spin-flip dynamics. The spin dynamics of cobalt assumes an intermediate position and changes weakly 
upon the transition from the bulk to 1 ML Co/Cu(100). 

Based on the examples above we suggest the following classification of the spin- wave Landau damping in 3d magnets. 
We distinguish three regimes: (A) for small wave vectors the acoustic magnon pole appears outside the energy range 
of the Stoner continuum and the spin- waves resemble closely atomic-spin like precession assumed in the Heisenberg 
model. For the elemental transition metal ferromagnets this region contains typically only a very small part of the 
Brillouin zone, whereas it might span large part of q space in the case of half-metals'^^; (B) in the energy region of low 
intensity Stoner excitations (Rcxks 3> ImxKS 7^ 0) the trend to form the coherent precession of atomic moments is 
still observed. However, it is opposed by the one-electron Stoner transitions at this energy. This hybridization results 
in a broadening of the spin- wave peak and, respectively, in finite life time of the magnon. The magnons, however, can 
still be regarded as well defined excitations. This regime is absent in the uniform electron model, but much of the 
spin-flip dynamics in Fe, Co and Ni falls into this category; (C) finally, "spin-wave disappearance" regime corresponds 
to the situation where the expected position of the magnon pole lies in the energy region with high density of Stoner 
states. Li this case RexKS ~ ImXKS Stiid no well defined resonance peak can form since the phase relationship between 
the external field and the induced exchange-correlation field is destroyed due to the intense excitations of incoherent 
Stoner pairs. 

In different systems the relative importance of different spin-wave damping regimes may vary. In this respect 
already the comparison of the bulk 3d metals brings interesting observations. In bcc Fe, along certain directions in 
the Brillouin zone, "spin-wave disappearance" regime sets in very abruptly. In fee Ni and fee and hep Co regimes A 
and B dominate and practically no spin-wave disappearance is seen. The case of ultrathin films can differ strongly 
compared to the corresponding bulk system under the infiuence of the dimensionality and the presence of non-magnetic 
substrate. 
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Appendix A: Representation of the spatial dependency of susceptibilities 

The susceptibilities depend on two spatial arguments and are represented in a separable basis as follows 

X(x,x') = ^XAA"y3A(x)(^A'(x')*, (Al) 
AA' 

where the functions ip\ form a complete basis, xaa' is the susceptibility matrix in this basis. 

In this work, we mostly used the following function set, referred to as Y-Ch basis, to represent the susceptibilities. 
The angular dependence in an atomic (Voronoi) cell s is cast into real spherical harmonics (f ) whereas the radial 
dependence is represented by Chebyshev polynomials £[)^ (^) 



(A2) 



25 



Here, = x — s^, where Sg is the center of the Voronoi cell containing point x. 8s (x) is the shape function equal 
1 when X is inside cell s and otherwise. In the atomic sphere approximation, the Voronoi cells are substituted 
by atomic spheres with radii Rs- Ts = Ir^j and = r^/r^. ^[a,6](^) is an invertible function mapping interval [a,b] 
into interval [—1,1]. The additional multiplier in (jA2p improves the convergence properties of the basis and is 
convenient when solving the susceptibility Dyson equation. The composite index A = sfilm determines a supermatrix 
structure of the susceptibility. 

The Y-Ch basis offers a complete, accurate, and efhcient representation of the spatial dependencies of the suscepti- 
bilities. Compared to other approaches^^^^® no assumptions regarding the nature of orbitals responsible for magnetism 
are necessary: all orbitals are included on an equal footing. Also the full spatial dependence of the exchange-correlation 
kernel is taken into account. The number of necessary Chebyshev polynomials per site per spherical harmonic needed 
for accurate representation of the susceptibilities for systems considered in this work varies between 8 and 16. The 
basis functions are localized on atomic sites and, unlike the plane wave basis, can be used equally well for representing 
the spatial dependencies in the periodic solids, at surfaces and interfaces and in finite clusters of atoms. 

In the case when the cells are approximated with spheres, the basis functions are orthonormal regarding slm 
indices but not the Chebyshev index fi. Prior to eigenvalue analysis it is convenient to transform the susceptibility 
matrices into an orthonormal representation. We use the Lowdin transformationi2^ii2£ based on the matrix square 
root algorithm of Denman and Beaversi^l. 

For systems featuring a discrete translational invariance, it is convenient to express the quantities using the following 
mixed r-q representation. We define the lattice Fourier transformation 

/(r,q)=^/(r + R)e->i-^, (A3) 

R 

where the summation proceeds over the crystal lattice and r belongs to the Wigner-Seitz cell ftws of the crystal, q is 
a vector in the first Brillouin zone flez- The retarded susceptibility x*"' (r, r', q) relates the components of the external 
field and the induced density (the frequency argument has been suppressed) 

Sn^ir,q)=J2[ drV^ (r, r', q)S^ (r', q), (A4) 

j O ws 

where 

X(r,r',q)=^x(r + R,r')e-"1-^. (A5) 

R 

The (r, r') dependency above is given in Y-Ch basis. 

The inelastic scattering experiments probe the imaginary part of the Fourier transformed susceptibility^'*, obtained 
by projecting it on the plane waves 

where r g f^vvs a-i^d K is the reciprocal lattice vector. The Fourier transformation is defined now as 

XKK'(q) ^ f^ws / drdr'e-'('i+K)-'-x(r,r',q)e'('i+K')-"-'. (A7) 

K,K' are reciprocal lattice vectors and q G r^BZ- The definition is consistent in the sense that in uniform systems 
one obtains 

XKK'(q)=x(q + K)5KK', x{q) = J d^x(^)e-'^- ^ . (A8) 



Appendix B: Evaluation of the integral in Eq. (I36|) and analytic continuation 

Our purpose is to transform the integral in Eq. p6p in such a way that in the complex integration contours as little 
as possible evaluations of KKRGF are performed close to the real axis. 

This can be easily achieved in the case when oj = and 7 = 2mrkBT = uj!^, n G as equation reduces to the 
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FIG. 19: Evaluating of the first (a) and the fourth term (b) in Eq. (|36p around the valence band. Crosses (x) denote 
the positions of the Fermionic poles 9n. Part of the integration is transformed into the sums over Fermionic frequencies 
(small circular contours around selected The integration weighted with the Fermi-Dirac distribution function^^^ can be 
approximated close to the chemical potential /i using Sommerfeld expansion. This way of integration allows to avoid evaluation 
of KKRGF close to its singularities originating from the valence states; the valence band is marked with thick orange line. The 
sum of the integrals a, b and the one given by Eq. (|B2|| gives joint contribution of the first and fourth term in Eq. (|36p . 



Matsubara susceptibility evaluated at Bosonic frequenciei 



,35.38 



(Bl) 



where 0„ = /i + iw^, ^ stands for the chemical potential and = (2n + l)7rfcBr is the fermionic Matsubara frequency. 
We remark that the temperature is introduced everywhere in this work for the computational expediency only, in 
order to smooth the discontinuity of the Fermi-Dirac distribution. T does not influence the results as long as it 
remains much smaller than the characteristic band width. 

The case a; is more difficult to handle. Here we use the procedure that we refer to as nearly real axis 
approach^^ We take advantage of the periodicity of friz) on the imaginary plane and consider only the case when 
,Z 9 M > 0. The ffi'st and the last terms of Eq. ([551) can be computed rather straightforwardly, since the 
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Green's functions involved are evaluated on the same complex semi-planes (upper in the case of the first term and 
the lower for the fourth term). The integration contours can be deformed in order to perform the integration away 
from the singularities corresponding to valence states, as presented in Fig. [191 For energies e well below the bottom 
of the valence band and above the core states Ec the respective contours can be taken parallel to the real axis; in 
this case their joint contribution reads 
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(B2) 



where A > |w/2|. In the vicinity of the core states the integration contours must be deformed again so the KKRGF 
evaluation is performed away from the singularities. The contribution given by Eq. (|B2p vanishes rigorously for uj = 
and 7 



In the case of finite lo it can still be safely neglected, providing that U-qT and w are small. 



The second and the third terms cancel each other only for a; = and 7 = 0;^. In the w ^ case they are quite 
intricate to compute, since they involve Green's functions evaluated simultaneously on both complex semi-planes. To 
minimize the computational effort we rewrite the joint contribution of the second and third terms as follows 
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(B3) 



The structure of the integral is presented in Fig. [20l it has been "symmetrized" with respect to the real axis. The 
evaluations of both GFs are performed now as far as possible from the Kohn-Sham poles. When resorting to Som- 
merfeld expansion, the integration is reduced to a finite range e € [~ f > ^] • The two Matsubara sums are computed 
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FIG. 20: The structure of the integral in Eq. (fB3ll . 




FIG. 21: The analytic structure of transverse susceptibility and schematic presentations of different analytic continuation 
schemes. For particular q the singularities of KS susceptibility form Stoner continuum (SC). The additional singularities 
introduced into the enhanced susceptibility by the Dyson equation, i.e. spin-waves, can appear outside the continuum; such 
magnons cannot decay via Landau mechanism. On the contrary, when the SW pole appears in the Stoner continuum it acquires 
a finite life time manifested by an apparent shift of the pole into the lower complex semi-plane. The analytic continuation in 
the case of the temperature susceptibility (a) involves the reconstruction of the real time dynamics based on the values given 
in points located on the imaginary axis; it is in general unstable, (b) Nearly real axis calculations lead to much better stability 
and accuracy of the analytic continuation. 



directly term by term. The integration involves now evaluating the Green's function at distance 7/2 from the real 
axis and for sufficiently large M rapid convergence can be achieved. 

It is well known that the analytic continuation poses a tricky and potentially unstable numerical problem;^^^ In the 
context of this work one faces in general two contradicting requirements. The evaluation of temperature susceptibility 
given by eq. (jBip . corresponding to the purely imaginary frequency in eq. p6p . is much easier to implement and 
numerically much faster. Unfortunately, the subsequent analytic continuation becomes pathologically unstable, since 
the distance between the points of the complex plane where the susceptibility is actually evaluated and the points on 
the real axis where we want to determine it by means of the analytical continuation is large, cf. Fig. 121b .. In fact, the 
continuation cannot be performed without assuming certain explicit analytical form of the susceptibility, as e.g. in 
Refi^. It is much more beneficial to work with oj ^ 0, Fig. 121b. In this case, the computational effort increases with 
decreasing 7, since much denser sampling in the Brillouin zone integration is necessary. Smaller 7, however, stabilizes 
greatly the subsequent analytic continuation procedure. 

To perform numerical analytic continuation, in most cases we employ a rational function (Fade) 
approximatioi>i22rJ^, where a complex function /(z) is represented by a ratio of two polynomials. Alternatively 
we resort to the method of Haas, Velicky and EhrenreichJ^iiS 



Appendix C: Explicit form of products of two KKR Green's functions 



By solving the multiple scattering problem in the atomic sphere approximation, the Korringa-Kohn-Rostoker (KKR) 
Green's function (GF), eq. ([55)1 . is obtained in the following representation 

G.(x,x',z) = yi^ij;i(r>,zXi(r<,z)<5™„ + ^<i(r™,z)G™£'^,(z)KL'(r;,z), (CI) 

L LL' 
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where cr = (a, /3) is the spin index. In the case of colhnear magnets, only the two diagonal spin components of 
the GF are non-zero (a = /? = fii)- r„ = x — s™, where Sm is the position of the atomic site closest to x, 
r< denotes the one of two vectors rm,rj„ being shorter and r> the longer one, i?™(r, z) = i?™(r, z)Yt(f) and 
H^{r,z) — R™{r)YL{r) are regular and irregular solutions of the radial Schrodinger equation for atomic site to. 
The first term in Eq. (jCl[) represents the single-scattering GF, while the second term describes multiple-scattering 
processes via the back-scattering operator G"^^,(z), which can be computed from the algebraic KKR Dyson equation 

G7El'{^^) = .9™'(^) + E 5E^"(^)i^L"(^)G^2"L' (C2) 

kX" 

where g™^, are the KKR structure constants and t™]^{z) is the single-site scattering matrix. We note that the 
computational method presented in this paper is trivially generalizable to the full potential treatment. 

In periodic systems, the product of two Green's functions appearing in Eq. (j36p involves additionally the convolution 
over the Brillouin zone 

S'^,^2(r,r',q, zi,Z2) = <flCrL / - — G^^ (r, r', k, zi)G^, (r', r, k - q, za). (03) 

D stands for the dimensionality of the periodic lattice. The number of necessary integration k— points decreases 
rapidly as one moves away from the singularities of KKRGF, i.e. Kohn-Sham energies located on the real axis. 
Upon substituting the KKR form of G in equation (jCSp we obtain 
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r"^G£;i,i3iJzi,Z2,q)<L,('^™,^l)<L.K,^l)^^.L3K,^2)i?r,L.('^n.,Z2) + 
'5™„'5L,L.V^<:L,(r<, zOi?;':L,(r>, zOSr.L3L4(^2)KU3(^™' ^2)i?r.L,(^m, ^2) + 

SrnnSL,LjL,L,V^2K,L, {r < , ZlW'^^L, ir> , Zl)i?™ ^3 (r< , Z,)H:1^^ (r> , Z2)). (C4) 

The first term comes from the convolution of two backscattering operators 

^"''crXL,LA^i,^2, q) = 75^ / d^kG;7^^^,(zi, k)G™^3^,(z2, k ~ q), (cs) 



while the next two terms involve only diagonal part of it 



B7lu{^^)^7^1 rf^kGXL'(^,k). (06) 



^BZ Jf2Bz 

By means of Gaunt coefficients the products of four spherical harmonics are reduced to pairs Yl{t)Yli{t'). The 
remaining radial dependence is approximated using Chebyshev polynomials. This gives the representation of the 
susceptibility in the Y-Oh basis. 
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